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SUMMARY 



This paper deals with the numerical solution of linear 
differential equations of fourth order by finite differences. 

It points out significant (but usually overlooked) errors which 
result from the conventional method of imposing the boundary 
conditions in such problems. Revised finite difference formulas 
are derived which apply near the boundaries and which eliminate 
the above errors. 

Three commonly encountered boundary conditions are consid- 
ered. These correspond, in the terminology of beam analysis, to 
a clamped end, to a simply supported end and to a free end. 

The improvement in accuracy that can be achieved with the 
revised formulas is illustrated by two representative examples. 
The revised formulas are shown to reduce the overall error of the 
numerical solution by a factor of about five in a typical case. 
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1. Introduction 



Let (x) denote an unknown function which is governed by 
a linear differential equation of the following form 

+ f^Cx) + f 2 (x) D^ 4 > + fj_(x) D(j) + f^(x) ^ = g(x) (1.1) 

where the coefficients f 2 (x), f^(x), fQ(x) and the 

forcing function g(x) are all known functions. 

We seek an approximate numerical solution of Eq. (1.1) by 
finite differences which satisfies appropriate boundary conditions 
at X = 0 and at x = i . Three commonly encountered bo'ondary 
conditions are considered in this paper, and are labelled below 
according to the terminology used in beam analysis, namely, 

1) Clamped End 



D4> 


-ol 


at 


o 

II 

X 


or 


at 


X = il 


(1.2) 


Simply 


Supported 


End 












<(> 


II II 


at 


X = 0 


or 


at 


X = 


(1.3) 


Free End 














2 

D (j) 


= 0, 


at 


X = 0 


or 


at 


X = i 


(1.4) 



3 

D (*) = 0 



The boundary conditions imposed at the two ends may be like 
or unlike. There is no restriction, except of course that if one 
end be free, the other is normally clamped to ensure that the 
beam configuration remains stable under arbitrary loading. 
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In the finite difference analysis, the domain 0<x < I 



is subdivided into N' equal intervals each of width 




If we include the two end points, the above domain contains 
(N' +1) stations. The function ()> is now represented by its 
discrete values (|)^ (i = 1, 2, 3, ....) at these stations. It is 

clear from Eqs . (1.2) and (1.3), however, that the fionction 

vanishes at either or at both of the end points depending on the 
particular boundary conditions that happen to apply. In any case 
it is convenient in the present context to let N denote the 
actual nximber of stations at which is initially unknown. For 

a beam whose ends are either clamped or simply supported, 

N = N' - 1. For a beam free at one end and clamped at the other, 

N = N' . 

It is also convenient to designate the stations at which <j> 
is initially unknown, from left to right, by index i = 1, 2, 3, 

. . . N . For a beam free at the left end, index i = 1 denotes 
the actual free end itself, at location x = 0 . For a beam 
clamped or simply supported at the left end, however, index i = 1 
denotes the first station inboard of the left end, at location 
X = h . This convention simplifies the indexing of the final 
matrix equations so that they always run from i = 1 to i = N 
inclusive . 

To determine the unknowns 'J* 3 ' • • • • rewrite 

Eq. (1.1) in finite difference form for the ith station, then 
require that the resulting equation be satisfied for each of the 
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stations i = 1, 2, 3, . . . ,N . This produces N simultaneous 

equations in N unknowns which suffices to establish the required 
solution. The actual manipulations are handled most expediently 
in matrix format. 

Except for the stations at or immediately adjacent to the 
boundaries, the various derivatives in Eq. (1.1) are usually 
approximated by the four conventional central difference formulas 
summarized in Eqs . (1) through (4) of Table I. 

These formulas are based on approximating the function ip 
in the vicinity of arbitrary station i by a truncated series of 
the form 

2 3 

2 (x-x. ) 2 (x-x.) 

(P - P-) = Dp. (x-x.) + D (J) . + D ij). + ... (1.6) 

The number of terms retained in this series depends on the accuracy 

required in the final difference formulas. Eq. (1.6) is applied 

to a number of contiguous stations symmetrically disposed on both 

sides of station i . This yields a set of simultaneous equations 

which can be solved for the initially unknown coefficients Dp^, 

2 2 

D D .... of the series thus yielding the required finite 

difference formulas. These formulas are represented by the 
bracketed expressions in Eqs. (1) through (4) of Table I. 

By re-introducing an additional term into the original series 
and retracing the above solution, an estimate of the truncation 
error is obtained. This is represented by the final term, that 
is, the term outside the brackets, in Eqs. (1) through (4) of Table 
I. Theoretically, the total truncation error could be represented 
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Table I Conventional Finite Difference Formulas 



Central Differences 



D0. = ij 



° - 5^i-l * “ ^ ^ 



1 2 3 

- i h^D-^0. + . . 

D 1 



( 1 ) 






0 + 0._3_ - 20. t 0.^^ + 0 



^ h^D"^0^ + . . 



( 2 ) 



D^0. = ^• 



i .0 . +0 + 0 “ 0 + —0 

2^i-2 ^i-1 ^ ^i+1 2- i+2 



- ihV0. + 

4 1 



(3) 






♦i-2 - “*i-l ^ + 2 



- i h'oS. 
0 1 



(4) 



Clamped End 



D0 = - 
^^1 h 



0 + ^2 + ^ 



i h^D^d, + 
6 1 



(5) 



D^0 = — 

1 h2 



- 20 ^ + 02+0 



h * 



( 6 ) 



D^0 = ^ 

1 



L 6 - 0 + i 0 

2 ^1 ^2 2 ^3 



°'*1 = ^ - '*2 * *3 



I °^*o * ■■■ 



1 



3 h 



(7) 



( 8 ) 
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Table I Cont'd 



Simply Supported End 



Ddi = ± 



0 + ^2 + 0 



g h D <^1 



D^di = ^ 



~ 20^ + (p 2 ■*■ 0 



T2 h D 01 + 






2^1 “ *^2 2^3 



- ^ hD%, + 



4 . 1 



+ 50^ - 402 + 0^ 



1 4 

^ h ^ 



Free End 



At station i = 1 






- ^1 + ^2 



- i +' . . . 

O i. 



D 0^ * 0 



D^0^ = 0 



d'^0 



1 w4 



20 ^ - 402 + 20 2 






At station i = 2 






■|0j^ + O+ -^3 + O 






^^2 1. 1 
° ^2 2 



01 - 20 2 + 03 + 0 



- ^ + 



^ 



0 + i^2 ■ '^3 + 



- k “'*1 ^ 



4 • 1 

D% = 



20 ^ + 502 “ 40 ^ + 0 ^ 



1 4 

+ ^ t 



(9) 

( 10 ) 
( 11 ) 
( 12 ) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

( 20 ) 
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by an infinite series of which the quantity shown is merely the 
leading term. 

Notice that all four of the conventional central difference 

2 

formulas shown in Table I have truncation errors of order h 

Additional terms of the series which represent the truncation 

error are not given in the table but they can be shown to involve 

3 4 5 

steadily ascending powers of h , that is, h , h , h , .... etc. 

Hence in the limit of very small mesh size h , the higher order 

become relatively negligible and the single term shown becomes 

itself an adequate approximation to the overall truncation error. 

2 

Also notice that for a truncation error of order h , the 

2 

central difference formulas for and D d> ^ require a band 

width of only three, that is, they involve the values of <f> at 

only the three successive stations (i-1) , i and (i+1) . To 

achieve the same order of truncation error for D'^d . and D^d • 

1 1 

on the other hand requires a band width of five. 

When reduced to matrix format, the basic relation Eq. (1,1) 

has a band width equal to the widest band of any derivative which 

appears in it. For a given order of truncation error, the deriva 

tive of highest order requires the greatest band width. Thus for 

2 

a truncation error of order h , the final calculation matrix 
representing Eq. (1.1) has a band width of five; it is a penta- 
diagonal matrix. 

The first four formulas in Table I may be applied routinely 
at interior stations, but special problems arise near the bounda- 
ries. It suffices here to illustrate the problem for the left 
end. Consider the case where the left end is either clamped or 
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simply supported. Then at station i=l immediately adjacent 
to the end in question, the conventional central difference 
expressions for D<^ , and D <(> , involve the three quantities 

^ 1 ' ^2 * corresponding central difference ex- 

• 3 4 

pressions for D , and D 0 involve the five quantities 

<t >2 and 0^ . The quantities <f>^, <p^ and 0^ 
entail no problems. The quantity vanishes for either a 

clamped end or a simply supported end and therefore involves no 
difficulties. The problem that now arises, however, is that the 
quantity 0_^ lies outside the normal domain of integration 
o < X < , and is in fact undefined. Thus D^0 , and D^0 , 

cannot be evaluated from the usual central difference formulas 
until this difficulty be resolved. 

A similar problem arises also if the left end be free. Re- 
call from the labelling convention adopted earlier that for a 
free end, the actual end of the beam at location x = 0 is 
denoted by station i = 1 . Application of the conventional 
central difference formulas at this station involves the five 
quantities ^ ^2' ^2 * Application of the conven- 

tional central difference formulas at the next station i=2 
involves the five quantities 0^, 0^, 02 / '^4 • Once again, 

however, the quantities 0_^ and 0^ are initially undefined 
and the conventional equations cannot be applied until this 
double indeterminacy be resolved. 

It is apparent from considerations of symmetry that 
corresponding questions arise also at the right end of the span. 
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Fortunately, there exists a generally accepted convention 
for resolving these questions. The rationale of this conventional 
solution, along with a careful analysis of its limitations, is 
presented in the next section. The results are in some respects 
surprising . 
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2 . The Conventional Boundary Formulas and Their Limitations 

The conventional solution for the left end clamped can be 
found by making the following substitutions in Eg. (1) of Table 
I, namely, 

i = 0 

0 ^ = 0 ( 2 . 1 ) 

D*o = “ 

The resulting relation can then be solved for 0 ^ as follows 
0_3_ = 0^3_ - j h^D^0^ + ... (2.2) 

We can now evaluate Eqs. (1) through (4) of Table I at station 
i = 1 making use of the foregoing substitutions for 0 and 0_, 
In this way we obtain the conventional finite difference formulas 
for a clamped end as summarized in Eqs. (5) through (8) of the 
table . 

A similar process may be used for the simply supported end. 
Thus we make the following substitutions in Eq. (2) of Table I, 
namely, 

i = 0 

0=0 (2.3) 

o 

D^0 = 0 

o 

The resulting relation can then be solved for 0_^ as follows 
4 > ^ = ~ *^ + l (2.4) 
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We can again evaluate Eqs . (1) through (4) of Table I at 

station i = 1 by making use of the above substitutions for 0 ^ 
and . In this way we obtain the conventional finite difference 

formulas for a simply supported end as summarized by Eqs. (9) 
through ( 12 ) of the table. 

The conventional solution for the left end free can be found 
by making the following substitutions in Eqs. (2) and (3) of Table 
I, namely, 

i = 1 

= 0 (2.5) 

= 0 

The resulting pair of equations can then be solved simultane- 
ously for 0^ and 0^ . This gives 



( 20 ^ -02 + 0 ) + 


^ h'^D'^ 0 ^ + . . . 


( 2 . 6 ) 


(40^ - 402 + 03 ) 


+ h'^D'^ 0 ^ + . . . 


(2.7) 



We can now evaluate Eqs. (1) through (4) of Table I at station 
i = 1 in the usual way, by making use of Eqs. (2.6) and (2.7) to 
eliminate the quantities 0^ and 0^ wherever they occur. In this 
way we obtain the conventional finite difference formulas at station 
i = 1 for a free end as summarized in Eqs. (13) through (16) of 
the table. 

Notice that for a free end, special formulas are required not 
only at station i = 1 but also at station i = 2 . At the latter 
point the quantity 0_^ is no longer involved so that Eq. (2.7) is 
no longer needed. But Eq. (2.6) is still needed to eliminate 0 
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wherever it occurs. In this way we obtain the conventional finite 
difference formulas at station i = 2 for a free end as summarized 
in Eqs. (17) through (20) of Table I. 

It is usually assumed that the conventional boundary formulas 

all have truncation errors of order h because they are derived 

from central difference formulas having truncation errors of order 
2 

h . Unfortunately, the results listed in Table I show that this 
is not necessarily the case. Notice that if the truncation error 
terms be omitted from Eqs. (2.2), (2.4), (2.6) and (2.7) by over- 

sight, the conventional difference formulas will all seem to be of 
2 

order h . This point is easily overlooked and it probably 
accounts for the widespread misunderstanding concerning the true 
truncation errors of the conventional finite difference formulas. 

An important aim of this paper is to call attention to this common 
oversight. 

It is clear from Table I that some of the conventional formulas 
have truncation errors of order h , h° or even h ^ ! Such errors 
are excessive and are inconsistent with the order of error involved 
in the rest of the calculation matrix. They can be expected to 
increase unnecessarily the overall error of the numerical solution. 
Sample calculations presented later in this paper confirm that 
this is indeed the case. It is shown that revision of these faulty 
expressions greatly reduces the overall error of the final solu- 
tion. 

The method of making the needed revisions, and the results of 
these revisions, are summarized in the next section. 
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3 . Revised Finite Difference Equations of Consistent Truncation 



Error 



Fortunately, the method of developing revised finite difference 
formulas where needed is straight forward, if sometimes rather 
tedious. The function <p in the vicinity of the end is repre- 
sented by a truncated power series. The form of the series must 
be such as to satisfy the boundary conditions of interest as 
previously summarized in Eqs . (1.2), (1.3) or (1.4) 

Thus for a clamped end at x = 0 we set 



<p = D 4> 



fr ^ °'^o 



3 

X 

3! 



4 

D 0 



TT 



4- 






(3.1) 



For a simply supported end at x = 0 we set 

* = D*^x + fr * • • • 



(3.2) 



For a free end at x = 0 we set 

.45 

(0 -0^) = D0^ X + D^0^ Iy + D^0^ fr 

The number of terms retained in the series depends as usual 

on the order of the derivative to be estimated and the truncation 

error of the estimate. The truncated series expression is applied 

to a number of contiguous stations near the end of the beam. This 

yields a set of simultaneous equations which can be solved for the 

2 3 

initially unknown coefficients D0 , D b , D 0 , .... or 

■“ o o o 

T 3 

D0^, D‘"0^, D <t>^ .... in terms of the quantities 0^, 02/ 0^/ ... . 
Once these coefficients are obtained in this way, the various deri- 
vatives of 0 at any specified station may be evaluated analytically 
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from the basic series expression. The result is the finite 
difference formula of interest. 

By introducing an additional term into the original series 
and retracing the original solution in the usual way, the leading 
terra of the truncation error can readily be found. 

The above process has been carried out for all of the con- 
ventional finite difference formulas of Table I which involve 

excessive truncation errors. The revisions are always such as to 

2 

render all truncation errors consistently of order h . The re- 
sults of this revision are summarized in Table II. Notice that 

every formula in Table I whose error term is initially of order 
2 

h reappears in Table II without change. Only the faulty formulas 
have been revised. 

An interesting detail appears in Eq. (8) of Table II. The 

2 

truncation error term is seen to involve h but the numerical 
coefficient happens to vanish. Hence the true truncation error 

3 

degenerates to order h in this particular instance, which is 

exceptional. It might be supposed that this circumstance would 

permit a reduction of band width here but such is not the case for 

if the band width be reduced by one, the resulting truncation 

2 

error is found to be of order h not of order h . Hence Eq. 

(8) must be retained as written. 
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Table II Finite Difference Formulas of Consistent Second Order 



Truncation Error 



Central Differences 



D(f) . = =■ 
h 



0 - + 0 + + 0 



- i . + 



u 



i ^ [0 + - 2« . + + 0 

n 



■ h ^ •• 



r^3^ • 1 

D 0i = ^ 



- I 0^_2 + <Ai_i + 0 - + I <^i+2 



- j h^D^(p. + 

4 1 



4 • 1 

0., = ^ 



<p . ^ - 40. , + 60. - 40. , + 0. 



i-2 i-1 



i+1 i+2 



- I- h^D^0 . + . . . 

O 1 



clamped End 



D0 = — 
^^1 h 



0 + i <>2 + 0 



- i . 



i ij 



-20., + 0_ + 0 



- n ^ 



i i- 

1 3 



- 301 t 0 t i 03 



3 v,2-5. , 

20 ^ ° ^0^ 



r^4A • 1 

D0, 



1602 - 902 + ! ^3 - I ^ 



+ (0) h^D^0 . + 
^1 



( 1 ) 

( 2 ) 

(3) 

(4) 



(5) 

( 6 ) 

(7) * 

( 8 ) * 
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Table II 



Cent ' d 



Simply Supported End 



“^1 ■ E 



0 + f <.2 0 



1 >,2n3A ^ 

- ^ h D 0^ + 






-20^ + 02 + 0 



T2 1 + 






Jl 01 - 11 *^2 ^ 11 *^3 



9 2 5 

+ ^ h^D^0 + 

44 o 






32. 27, .8, 1, 

^ 01 - ^ 02 + 5 03 - To 04 



2 ^2^6^, + 
- ^ h D 0Q 



Free End 



At i = 1 



_ 1 
1 h 



D0T=n- |^“'^i‘^'^2’^^ 



1 2 3 

- ^ h D-^0, + 
6 1 



D 0^ = 0 



D 0^ = 0 



4 • 1 



4.23530^ - 9.176502 + 5.647102 - 0.70590^ 



+ 0.2255 h^D®0^ + 
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Table II Cont'd 



At i = 2 






1 3 

- I h“D^02 ^ 



D^4>. 



^ M<. - 



2^2 + + 0 



- ^ hV., . 



°'*2 =^ 1 - 15 *! - 1*3 + “ 



* h * 



D^4> = — 

2 h-' 



- 1.64710^ + 4.235302 " 3.529402 + 0.94120^ 



+ 1.1451 h2o^0^ + 



(17) 

(18) 

(19) * 

(20) * 



*The asterisks denote those formulas of Table II which differ from 
the corresponding formulas of Table I. 
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4 



Finite Difference Formulas of Minimum Truncation Error for 



Specified Band Width 

An interesting and significant feature of Table II is that 
the finite difference expressions for the various lower derivatives 
usually involve lesser band widths than that required for calcula- 

, 4 

ting D 4 > . The latter is fixed by the requirement that the 

. . 4 

truncation error in the approximation for D 0 should be of order 




A natural question that now arises is whether there be any 
possible further advantage in making use of this full available 
band width to reduce the truncation errors associated with these 
lower derivatives. 

This approach is attractive in that it entails no increase 
whatever in the overall band width of the final calculation matrix. 
Hence it can be accomplished with a negligible increase in the 
calculation burden. 

On the other hand, the method has the limitation that, although 

the accuracy of the lower derivatives is improved, the error in the 

2 

overall equation still remains of order h . Hence what we have 

here at best is the somewhat uncertain possibility of a marginal 

further improvement in accuracy at essentially negligible cost. 

In order to permit investigation of this possibility, a further 

revision of the finite difference formulas was therefore carried 

out based on two principles. The first was to choose a band width 

4 2 

such as to approximate D 0 to order h . The second was to 
utilize the resulting full band with for estimating all lower 
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derivatives. The finite difference formulas obtained in this way, 
along with their associated truncation errors, are sxammarized in 

4 

Table III. Notice that all expressions for D <p are the same in 
Table III as in Table II. Only the lower derivatives are different. 
A further significant feature of Tables II and III is that 

4 

the expressions for D <p ^ at station i=l involve terms in 
0^, 0^/ ^3/ ^4 • The first three of these fall within the overall 
penta-diagonal matrix format which characterizes all interior 
stations i = 2, 3, 4, ... etc. But the term in 0^ falls outside 

this band at i = 1 . 

Similar considerations also apply at the opposite end, station 
4 N 

1 = N . Thus D 0 invloves terms in ^^^ 2 ' ^N-2' "^N-l 

and the term in ^^^.3 falls outside the basic penta-diagonal 

format. 

Fortunately, there is a very simple method to eliminate the 
term in 0^ from the equation at i = 1 and the term in ^^-2 
from the equation at i = N . This restores the strictly penta- 
diagonal format which is so convenient and efficient in the final 
numerical solution. The details of the method are explained in a 
later section of this paper. 
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Table III Finite Difference Formulas of Minimum Truncation Error 
for Specified Band Width** 



Central Differences 



= K n *i-2 - I *i-i *1 



12 ^i +2 



+ ij * 






D^ 0 . = Kr 
^ h-^ 



I2 ^i-2 + 3 -^i-l 



5 a ^ 4 , 

2 0 i + 3 



1 A 

- 12 ^i +2 






j- I *i-2 • *i-l ^ ^ k *i«l - I * 



^ r *i-2 ■ “*i-l 



- i h^D^ 0 . + 
6 1 



Clamped End 



I "l * I '^2 



6 ‘*'3 ■*■ 43 *4 



no 



= n j- r ^ * I *2 - f ^3 * 5 

K *• 



48 



3 k 



n3A 1 1 



° ■ 7 *2 * i '^3 ■ 16 *4 






4 • 1 

0 



16*1 - 9*2 * 1 * 3 - 1*4 



+ ( 0 ) + 



(D* 

( 2 ) * 

(3) 

(4) 

(5) * 

( 6 ) * 

(7) * 

( 8 ) 
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Table III Cont'd 



Simply Supported End 



D01 - ^ 



47 . ^93 , 

150 ^1 100 ^2 



11 ^ J. 17 A 
50 ^3 600 ^4 



9000 ^ ^ 



D^0 = i_ 



38. ,29, 2..1, 

T5 "^l 20 "^2 ■ I5 "^3 I20 *^4 



^ mo * 



„3, ^ 1 12 , 63 . , 8 , 9 . 

° - 17 [— ■ 20 '^2 ■" 5 *3 ■ 40 *4 



, 558 v3„6, , 

■" 5400 h D + ... 



-4, j, 1 f32 , 27 , 

° *1 - 74 ( 5 - '^1 ' ?- * 



8 



2 5 ^3 10 "^4 



2 ^2^6 



C ^ “ 1 A 71 I “ AC ^ ■*■ 



25 



Free End 



At station i = 1 



= E 



- 1.12740^ + 1.270602 



0.158802 + 0 . 01570^1 

- 0.0035 h^D^0j_ + 



D201 = ° 



D-^0^ = 0 



r^4A • 1 



4.23530^ - 9.176502 + 5.64710^ “ 0.70590^ 



+ 2255 h^D^0^+ 



At station i = 2 



^2 = h 



— (h + + — 0 — i— ^ 

3 ’^l 10 10 3 30 4 



+ 0.0106 



Dh = 



1.1373<^^ - 2.3529<^2 + 1.29il<t> ^ - 0.0784<^^| 

+ 0.0740 + 
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Table III Cont'd 





1.29410^ 



2.470602 1.058802 + 0,11770^ 

+ 0,3765 h^D^0^ + . . . 




1.64710^ + 4.235302 " 3.529402 + 0.94120^ 

+ 1.1451 h^D®0^ + . . . 



*The asterisks denote those formulas of Table III which differ 
from the corresponding formulas of Table II. 

**In this table the band width is always chosen so as to give 

2 4 

a truncation error of order h for D 0 . The truncation errors 
of the other derivatives are then fixed accordingly. 



( 19 ) 



( 20 ) 
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5. Difference Formulas at Right Boundary 



It has been convenient to present the boundary formulas in 
Tables I, II and III specifically for the left end of the span. 

The corresponding formulas for the right end can readily be in- 
ferred from the foregoing results from considerations of symmetry. 
First we replace 4 > 2 > ^ 3 ^ ^rid by the corresponding quanti- 
ties 0^, "^N-1' "^N-2 '^N-4 -^sspectively . In this connection 

it is also convenient to reverse the order of corresponding terms 
from left to right so that the terms still appear in order of in- 
creasing index. Secondly, we must reverse all signs for the 
derivatives of odd order. 

We illustrate these rules by applying them to a clamped end, 
as given originally in Eqs . (5) through ( 8 ) of Table I. Thus for 

the right end clamped 





(5.1) 




(5.2) 




(5.3) 




(5.4) 



The same method can of course be applied to any of the other 



results developed in this paper to adapt them to the right end of 
the span . 



6 . Reduction to Penta-Diagonal Format at Boundaries 



If the basic relation, Eq. (1.1), be expressed in terms of 
the conventional finite difference formulas of Table I, the re- 
sulting calculation matrix is strictly penta-diagonal . It has 
earlier been pointed out, however, that if Eq. (1.1) be expressed 
in terms of the revised finite difference formulas either of Table 
II or of Table III, the resulting calculation matrix, while still 
penta-diagonal over all interior rows, falls outside the penta- 
diagonal format by one place in the first row i = 1 and last row 
i = N . The purpose of this section is to show how these two 
minor deviations can easily be rectified. 

Let us consider the left end of the span first. More speci- 
fically, consider the final matrix equations at stations i = 1 
and i = 2 . These will reduce to the form 



^11 *^1 


“12 


4>2 + 


®13 


(j) 3 + 


“14 *4 


f-i 

ii 


(6.1) 


“21 ^'1 ^ 


“22 


4>2 + 


“23 


^3 + 


“24 *4 


^2 


(6.2) 



where the coefficients and the quantities g^ and 

are all known constants in any particular case. 

We can eliminate 9 ^ between these two equations by multi- 
plying Eq. (6.1) through by «24 (6.2) by - , then 

adding. These operations give rise to the following auxiliary 
constants which can be readily computed as definite numbers, namely, 

“11 " “24“ll " “14®21 
“12 = “24“12 ■ “14“22 
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“l3 “24“l3 “ “l4“23 

I 

~ ^24^1 ~ °*14*^2 (6.3) 

In this way we finally obtain the relation 

III I 

ot^2^<^2^H“ct^2*^2"^^13^3 (6.4) 

Eq. (6.1) of the original matrix can now be replaced by Eq. 
(6.4) . This substitution reduces the band width at station i = 1 
such as to fit the matrix into the required penta-diagonal format. 

A corresponding correction can also be made at station 
i = N . 

Notice that these operations , which reduce the band width in 
the first and last lines of the final calculation matrix as re- 
quired, are very simple and do not adversely affect the truncation 
error of the final system of equations. 
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7. Example; Uniformly Loaded Beam with Clamped Ends 



Two examples have been worked out to illustrate the improve- 
ment in accuracy that occurs when the revised difference formulas 
of Table II or Table III are used in place of the conventional 
formulas of Table I. In this section we consider a beam of uniform 
stiffness with clamped ends under a uniform distributed load. In 
the next section we consider a beam of variable stiffness with 
simply supported ends under a uniform distributed load. 

For the uniform beam, the nondimens ional governing equation 
reduces to 

= 1 (7.1) 



and the boundary conditions are for clamped ends 



0 = 0 
D0 = 0 



at X = 0 and x = 1 



(7.2) 



The exact analytical solution of these equations is simply 



^ (l-x)2 



(7.3) 



In the finite difference solution, only half the beam need be 
considered because of symmetry. Let index M denote the station 
at mid-span. The central difference expressions at stations 
i = M-1 and i = M may be simplified by the substitutions 



"^M+l ^^M-l 






(7.4) 
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The conventional finite difference matrix for this case can 
now be written almost by inspection. For the first row i = 1 
we simply utilize Eq. (8) of Table I. For the interior rows 
i = 2, 3, 4, ... we use Eq. (4) of Table I. The matrix coeff- 
cients for the last two rows i = M-1 and i = M are easily 
determined by making use of Eqs. (7.4) in connection with Eq. (4) 
of the table. 

The following matrix is thereby obtained. 




(7.5) 



The matrix for the revised difference scheme of either 
Table II or Table III is identical to Eq. (7.5) except for the 
first row. Using the revised difference formula shown in Eq. (8) 
of Table II or Table III gives the following improved approximation 
for the first row, namely, 



^ [1601 - 5^2 •^1^3 - 1*4]= 1 



(7,6) 
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The second row remains unchanged, namely, 




40 ^ + 602 ■ 40 ^ + 0 ^ 



1 



(7.7) 



Multiplying 
(7.6) eliminates 



Eq.(7.7) through by + j 
0. and gives 



7 



150 - H 0 + 1 0 

^^1 2 ^2 3 ^3 



and adding it to Eq. 




(7.8) 



Eq. (7.8) is the contracted version of the revised first row 
that is inserted into the final matrix. This version is just as 
accurate as Eq. (7.6) but has the advantage that it fits strictly 
within the penta-diagonal format. 

Solutions have been obtained in double precision to Eq. (7.5) 
for both the conventional first row and for the revised first row 
given by Eq. (7.8). Results were calculated for three values of 
the important dimensionless mesh size parameter h , namely, 
h = 1/10, h = 1/20 and h = 1/40 . The calculations were per- 
formed in double precision on the IBM 360/67 computer at the Naval 
Postgraduate School. 

As convenient measures of the overall relative error of the 
final numerical solution, the following two parameters are used, 
namely. 



and 



e 




max 




x=l/2 



2 2 
D d> - D ( 



6 = 



max 



x=0 



(7.9) 



(7.10) 
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In beam problems £ expresses the maximxmi relative error 
in the calculated deflection while 5 expresses the maximum 
relative error in the calculated curvature, and hence in the 
associated bending stresses. 

Curves of e and 5 versus h are shown in Fig. 1 corres- 
ponding to the conventional difference formulas of Table I and to 
the revised formulas of Tables II and III. 

The results are quite interesting. They confirm that both 
error indices e and 6 are greatest for the conventional difference 
formulas of Table I. At the other extreme, both e and 5 vanish 
identically for the difference formulas of Table III. Of course, 
such a complete elimination of all truncation error must be re- 
garded as exceptional. It can only happen for those cases, like 
the present example, in which the true deflection curve happens to 
be a quartic. The curve is a quartic for any uniform beam under 
uniform load, irrespective of the boundary conditions. 

The revised finite difference formulas of Table II are seen 
to produce results in this case intermediate between those given 
by Tables I and III. Thus e vanishes but 5 does not. By an 
interesting coincidence, the formulas of Tables I and II happen 
to give identical curves of 6 versus h in this example. 

Inspection of the detailed solutions reveal, however, that the 
actual truncation errors represented by Table II are precisely 
equal in magnitude but opposite in sign to those produced by the 
conventional formulas of Table I! No specific reason is known for 
this curious and striking co-incidence. 
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RELATIVE ERRORS e AND S 




h — > 



FIG. ERRORS OF NUMERICAL SOLUTION 
FOR BEAM OF UNIFORM STIFFNESS 
WITH CLAMPED ENDS. 



8« Example: Beam of Variable Stiffness with Simply Supported Ends 

In this section we consider a uniformly loaded beam of vari- 
able stiffness, simply supported at both ends. Specifically, we 
choose a beam governed by the nondimensional equation 

232 34 23 2 

D (xD^.)=xD(}) + 6x D 4) + 6xD 4) = 1 (8.1) 

over the range 

1 < X < 2 (8.2) 

The boundary conditions for simply supported ends are 

(j, = 0) 

> at X = 1 and x = 2 (8.3) 

D^<j) = o' 



The exact solution for 4> is 



j (10 Jin 2-3) (1-x) + i 



— + ( 3+x) £ 

X 



n X -xj 



(8.4) 



Next we must obtain the conventional finite difference 
matrix corresponding to Table I. It is convenient in this case to 

4 

multiply Eq. (8.1) through by h and to introduce the notation 

x^ = 1 + i h i = 1, 2, 3, .... (8.5) 

For all stations other than station i = 1 and station 
i = N, the regular central difference formulas of Eqs . (2), (3) 

and (4) of Table I apply. Upon substituting these equations into 
Eq. (8.1) and regrouping terms, we can reduce the result to the 
following format, namely. 
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(xj - 3x?h) <j>- . + (- 4x? + 6x?h + 6x.h^) <t>. , 

-L X X X X**X 

+ (6x^ - 12x.h^) 0. + (- 4x^ - 6x?h + 6x.h^) 0 . ^ (8.6) 

X X X X X X x + 1 

+ (xj + 3xjh) 0.^2 = 

For station i = 1 , Eqs . (10), (11) and (12) of Table I apply 

and Eq . (8.1) reduces in this case to 

(5x^ + 3x^h - 12x^h^) 0^ + (- 4x^ - 6x^h + 6x^h^) 4> ^ (8.7) 

+ (x^ + 3x^h) 0^ = h^ 

2 

For station i = N , expressions must be obtained for D 0^^ , 
D^0j^, D^0j^ by analogy with Eqs. (10) , (11) and (12) of Table I, 

using the rules of symmetry summarized in section 5. Then Eq. 

(8.1) reduces to 

(Xn - 3Xyh) + (- 4x^ + 6x^h + (8.8) 

+ (5x^ - 3x^h - 12x^h^) 

The finite difference matrix defined by Eqs. (8.6), (8.7) 

and (8.8) was solved nxamerically for various values of the dimen- 
sionless mesh size parameter h . Each result was compared with 
the corresponding exact analytical solution of Eq. (8.4) for the 
same value of h . The relative error e , as defined by Eq. (7.9) , 
was calculated for each solution. The second error function 5 , 
as defined in Eq. (7.10), was not calculated for the present 
example . 

The above values of e were plotted against h . The re- 
sult is represented in Fig. 2 by the line marked I. 
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Consider next the revised solution based upon the revised 
finite difference formulas proposed in Table II. Careful examina- 
tion of Eq. (8.1) and of Tables I and II discloses that this case 
differs from the previous one only in the first and last rows, 
i = 1 and i = N . 

The revised equations for the first and last rows turn out 
to be as follows. 



(|^ xj + ^ x^h - 12 x^h^) 

+ (- xj - xjh + 6x^h^) 4>2 
+ (f xj 4. 30 03 + ^ x3) 0^ = h^ 



(8.9) 



and 



10 "^N-3 ^ ^5 ■ 11 "^N-; 



27 



54 



32 



IT ^ ’^N-l ^5" 






ii x^h 
11 V 



( 8 . 10 ) 

12x^h2)0^ = h^ 



The resulting matrix was then completed, solved and analyzed 
in the same general fashion as described above for the previous 
case. The corresponding curve cf error e versus mesh size h 
is shown in Fig. 2 by the line marked II. 

The same general pattern was also repeated for the finite 
difference formulas of Table III. In this case the entire matrix 
must be modified. The new matrix is defined by the following 
three governing equations. 
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RELATIVE ERROR 




FIG. Z- ERRORS ^OF VARIOUS APPROXIMATIONS 
FOR BEAM OF VARIABLE STIFFNESS 
WITH SIMPLY SUPPORTED ENDS. 
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For all interior rows 2<i ^ (N-1) the basic equation is 
found by substituting Eqs . (2), (3) and (4) of Table III into 
Eq. (8.1) This reduces to the result 

(xj - 3xjh - i x.h^) 0. , + {- 4x? + 6x^h + 8x.h^) 0. , 

1 12 1 1-2 1 1 1 ' 1-1 

+ (6x? - ISx.h^) 0. + (- 4x? - 6x?h + 8x.h^) 0. , 

1 11 11 1 1+1 



,.3^-2, 1 ,2, , , ^ 

+ (X. + 3x.b - 2 ^ih ) ^i+2 = h 



( 8 . 11 ) 



The result for row i = 1 is found by substituting Eqs. (10), 

(11) and (12) of Table III into Eq. (8.1) . This gives 

.32 ^3 , 72 2 , 76 „ . 2. , ^ , 27 3 189 2 , ^ 87 .2, , 

5 ^1 5 ^1^ " 5 ^1^ ^ ^1 ■ 5 ^1 “ 10 ^1 ^ 10 ^1^ ^ *^2 



^ ^8 ^3 ^ 48 2, 4 u2, ^ / 1 3 

+ (_ + 3- x^h - 5 x^h ) <^>3 + (- lo ^1 



27 2. ^ 1 ,2, , 

20 ■" 20 ' *4 



= h 



( 8 . 12 ) 



The corresponding result for station i = N is found from the 
above by applying the rules of symmetry as explained in section 5 . 
This gives 



10 ^ 20 ''' 20 ^ “^N-3 ^5 ^ " 5 ~ 5 ^ "^N-: 



+ (- 



27 



189 



10 



87 



32 



10 ^ ^N-1 ''' ^5 



72 ^2 76 2 

5 N 5 ’ 



N 



= h 



(8.13) 



Of course, the above results for rows i = 1 and i = N 
now iie outside the band limits for a penta— diagonal matrix, but 
they can readily be contracted by the method explained in section 

6 . 
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Eqs. (8.11), (8.12) and (8.13) were solved momerically in 

the usual way. The resulting curve of e versus h is shown 

in Fig. 2 as the line marked III. The marked change in the error 

curve between h = 0.1 and 0.05 is due to the fact that the 

sign of the error changes at the centerline. 

An interesting comparison of the results from the methods 

proposed herein can be made with the results of a study using the 

"half-station" method presented in Ref. 1. In the half-station 

method the finite difference approximations are made before ex- 

2 3 2 

pending the derivatives in D (x D 0) as opposed to the conven- 
tional, or whole-station method, wherein the finite difference 
. 2 3 2 

approximations are made after D (x D 0) is totally expanded. 

The conventional finite difference approximation for the boundary 
2 

condition D 0 was used at the two ends of the beam in Ref. 1. 
The relative error e , as defined by Eq. (7.9), produced by the 
half-station method for this example was shown in Ref. 1 to be of 
the form 

e = (0.4 X 10"^/0.4196 x 10~^) h^ (8.14) 

-2 

where the value 0.4 x 10 is obtained from Fig. 3 of Ref. 1 

-2 

and 0.4196 x 10 is the displacement of the beam at the center- 
line. This expression for e was verified by a niimerical 
solution of the half-station equations for several values of h . 
The error given by Eq. (8.14) is plotted in Fig. 2 as the line 
marked H . 
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Fig. 2 permits some very interesting and useful conclusions 
to be drawn. Like the previous example, it again confirms that 
the conventional finite difference formulas of Table I give by 
far the highest overall errors. The revised formulas of Table 
II are seen to reduce the above errors by a factor of approxi- 
mately five. It is interesting that curves II and III cross at 
approximately h a 0.05 , and that for mesh sizes finer than 
this, the formulas of Table II give somewhat better results than 
those of Table III. Inasmuch as the former are also somewhat 
simpler, there would seldom appear to be much reason to resort 
to the latter. 

Notice too that while the half station method of Ref. 1 is 
superior to the conventional formulas of Table I, the revised 
formulas of Table II are in turn superior to the half station 
method. 

We may conclude that the revised finite difference formulas 
of Table II probably represent the best overall compromise between 
accuracy and simplicity for most typical problems. 
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